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Abstract. We calculate cosmic ray transport with a collisional Boltzmann Transport 
. Equation, including E-M forces. We apply the Crank-Nicholson Method to solve this equa- 

HH ' tion. At each time step, the spatial distribution of cosmic rays is applied to the ZEUS 2D 

' magnetohydrodynamics model, which is then utilized to calculate the resulting E-M field. 

^ , Finally, the field is applied to the initial equation. This sequence is repeated over many time 

CL(' steps until a steady state is reached. We consider results from t = until steady state for an 

isotropic low energy cosmic ray flux, and also for an enhanced cosmic ray flux impinging on 
one side of a molecular cloud. This cosmic ray flux is used to determine an ionization rate of 
interstellar hydrogen by cosmic rays, f . Astrochemical implications are briefly mentioned. 
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^ ■ 1. Introduction slope of ~ E (Spitzer & Tomasko 1968, pre- 

' dieting very few low energy cosmic rays). A 

^ ■ Low energy (< 1 GeV) cosmic rays drive inter- ^^^^^^ ^^^^^ ^^^^-^ ^^^^^ 

stellar chemistry and may cause specific spec- ^^^^^^ necessary in order to under- 
tral features recently measured, such as the ^^^^^ ^^.^ flux-spectrum as a function of posi- 
6.7 keV emission line. Yet the origin and flux ^-^^ ^-^^-^ ^ molecular cloud. Modelling low 
of low energy cosmic rays is currently un- ^^^^-^ streaming will afford bet- 
known because the Sun s magnetic field de- understanding of interstellar chemistry and 
fleets these particles, so that they cannot be di- ^^j^ j^^^^ emissions caused by these cosmic 
^ ■ rectly observed. There is a great deal of un- ^^^^^^ ^^^^ ^ ^^^^^ ^^^^ 
rS _ certainty about the correct cosmic -ray flux- 
^ . spectrum for low energy cosmic rays, rang- 
ing from a steep slope of ~ £^ (predicting a 2. Equation and Numerical Method 

great many low energy cosmic rays, see Nath 

& Biermann (1994) all the way to a positive We have absolute coordinates xq and yo defin- 

ing the position within the cloud. For each 
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pointed in directions determined by angle 
we define a distribution function for cosmic 
rays, /. In the scheme utilized by Skilling 
(1975), we solve / using a semi-collisional rel- 
ativistic Boltzmann Equation in a two-fluid ap- 
proximation, where the cosmic rays are treated 
as one fluid, and the interstellar medium as an- 
other fluid. The equation is semi-collisional in 
that collisions within the interstellar medium 
are treated, as are collisions between cosmic 
rays and the medium, but not cosmic rays col- 
liding with other cosmic rays. 

In our two dimensional scenario, we set up 
a local coordinate system (x,y) for each small 
area, such that x • B = B and y • B = 0. B 
is separated from xo by an angle a such that 
the simple rotational transformation will map 
xo.yo x,y. We then apply the Fokker-Planck 
equation to solve for the distribution function. 
We use the Fokker Planck equation of a form 
similar to that of Cesarsky & Volk (1978). 

df c/xp df \ - jj} cp d\a^B df 
dt y dx 2 y dx 5yu 



dp 



/.^(x.v)§^.i^(y.v)f 

ox 1 oy 



df I ^ df 
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ox 2 ay 



^( dXoudi ^ di^i^" di) 2(dp) ' 
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Now, dfldt = dfldt + (v • V)/, c the speed of 
light, y - (I - v^/c^y^^^, is the pitch angle 
diffiision coefficient, and u is the velocity of a 
mean wave frame, or v -n va. The l.h.s. terms in 
Equation (1) describe the change in motion of 
relativistic charged particles due to convective 
motion of the plasma itself, as well as cosmic 
ray streaming along the magnetic field. The 
first two terms on the r.h.s. describe the eff'ects 
of the change in plasma densities. The third 
term is momentum change due to inelastic and 
elastic two-body collisions with the medium, 
and the fourth term describes the scattering of 
particles due to two-body collisions as well as 
magnetic field irregularities. The last term on 
the r.h.s., and the largest change we have made 
to Cesarsky & Volk (1978), is an addition from 



SkilUng (1975). This term accounts for spatial 
difliision across B. 

This equation does not include source 
terms or acceleration mechanisms, although 
they are relatively straight-forward to include 
both in Equation (1) and its numerical solution, 
so long as some already-estabUshed accelera- 
tion mechanism is provided. 

For collisions, we separate the coUisional 
momentum change into elastic and inelastic 
terms, referred by the subscripts "in" and "el", 
respectively, and the approximation for the in- 
elastic case is: 



(dp\ 
[ dt /„ 



ym 



Ap; 



(2) 



where o-, is the inelastic scattering cross- 
section from Cravens et al. (1975) and other 
sources, listed and reviewed very well in 
Padovani et al. (2009). The other terms, Ap 
is the momentum change from each collision, 
also reviewed in Padovani et al. (2009) and 
Rimmer et al. (2011). n is the density of the 
cloud, and m is the mass of the cosmic ray par- 
ticle, either the electron or proton mass. Elastic 
scattering is dealt with in a similar manner, ex- 
cept that the momentum is conserved over the 
two bodies involved in the collision, and the 
scattering cross-section is difl'erent. It is impor- 
tant to note that the elastic scattering also im- 
pacts D^. 

We solve this equation using the 
Crank-Nicolson method (Crank et al. 
1947), evolving the system from 
xi,yi,pi,iii,ti -> Xi,yi,pi,fii,ti+i, and aU 
iterations thereof. We approximate the first 
and second derivatives to, for example: 

d£ ^ fjxj, yi, Pi, fij, tj+i)- fjxj, yi, pi, m, tj) 
dt ~ 



At 



d^f 1 



- 2(Ax)4^^^''^^'^^'^'-^"^'-^'^ 
-2f(xi,yi,pi,Hi,ti+i) + fixi,yi,pi,Hi-uti+i)j 

+{fixi, yi, Pi,fii+i , ti) - 2f{xi, yi. Pi, Hi, ti) 

+fixi,yi,Pi,i^i-i,ti)j; 

(3) 
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where A? and Ax are characteristic time and 
length scales. For a cloud one parsec in di- 
ameter, with consideration for the constants 
in Equation (1), the At must be set to less 
than ~ 1 year without becoming too inaccu- 
rate. What is advantageous about the Crank- 
Nicolson method is its stability. The method 
will not lose stability pretty-much regardless of 
the choice of length and time scales, and so it 
is relatively straight-forward to determine the 
self-consistent accuracy of the calculations. 

All the values for / for x,y,p,^ ai t - 
are given, as are the values for / at the bound- 
aries, X = Q,L;y = Q,L over all values of p,n. 
In this proceeding, we use a set flux-spectrum 
to determine the value for / at the boundaries, 
and this value does not change with time. At 
; = 0, the value for / is determined by the flux- 
spectrum at the boundary, and / = at aU other 
points in the cloud. 

The calculation proceeds from the initial 
conditions of Xi,yi,pi,fii,to = 0, applying the 
values in Equations (3) to Equation (1), ar- 
ranging the elements as a series of matrices, 
inverting these matrices, and solving for the 
unknown values of / at x,,y,, p,,//,, ti. This is 
continued until steady-staet is reached. 

At the same time, the electromagnetic field 
and local density are determined using the 
ZEUS magnetohydrodynamics code (?). The 
input to the ZEUS code is a charge distribution 
provided by the cosmic ray distribution func- 
tions for electrons and protons. The results of 
the ZEUS code are applied repeatedly for each 
time step for the transport equation. 

3. Results in terms of the Ionization 
Rate 

It is useful for astrochemists, and also con- 
ceptually advantageous, to represent the two- 
dimensional results for the cosmic ray dis- 
tribution in terms of a cosmic ray ionization 
rate, f which is the rate at which hydrogen 
atoms are ionized by cosmic rays. This can 
be achieved mathematically by converting the 
distribution function to a position-dependent 
flux-density, j{x,y,p) - p^f(x,y,p). To derive 
a position-dependent ionization rate from the 
flux-density, we use the form from Spitzer & 



Tomasko (1968) with constants in front to ac- 
count for ionization caused by the products of 
the first ionization, xi, discussed in Dalgarno 
et al. (1999). This is the ionization rate for pro- 
tons: 

j(x,y,E)(ri,pdE. (4) 

In this equation, cr, is the ionizing cross-section 
and iiinin is the minimum energy for cosmic ray 
ionization. We achieve f by averaging f along 
a fine passing through the two-dimensional 
cloud. 

We performed the calculations for 
f(x,y, p,ij). For both calculations, the flux at 
the boundary is taken from Nath & Biermann 
(1994). For the first case, the flux is isotropic 
and there is a low-energy cutoff for the initial 
flux-density of 1 MeV (of course, the flux 
density inside the cloud can extend down 
to fimin)- In the second case, the initial flux 
extends down to Emin but impinges only on 
one side. The other side has the same initial 
flux-density, but with the 1 MeV cutoff. 

For the first case, the cosmic ray ionization 
rate extends from about 7.5 x 10~" s~' at the 
center to 2 x 10""' s"' at the edges. This dif- 
ference is too small to accurately detect, given 
that chemical tracers are the best current way 
to determine the cosmic ray ionization rate, 
and are accurate only to within a factor of 2 
or 3 (see McCafl et al. 2003; Indriolo et al. 
2007; Le Petit et al. 2006, for a review). In 
the second case, however, the ionization rate 
spans two orders of magnitude, and should def- 
initely be within detection capability, provided 
that sources can be found near the sites of cos- 
mic ray production and with angular resolution 
capable of achieving length-scales of about 10- 
100 AU. 

4. Discussion and Future Woric 

To more thoroughly examine the ionization of 
cosmic rays, we need to treat electrons as well 
as protons. The cross-sections have already 
been included in the code, and the electron cos- 
mic ray streaming will be calculated simulta- 
neously with the proton cosmic rays as a logi- 
cal next step. EventuaUy a third dimension and 
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Fig. 1. ^ as a function of depth into the cloud 
for an isotropic flux at the boundary from Nath & 
Biermann (1994) with a minimum energy of 1 MeV. 
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Fig. 2. ^ as a function of depth into the cloud 
for an isotropic flux at the boundary from Nath & 
Biermann (1994) with a minimum energy of 1 MeV 
on the right side, and 100 eV on the left side. 



turbulence, as well as self-gravitation will be 
incorporated in the calculations. 

There are many other questions such a 
model may answer beyond the cosmic ray 
ionization rate, such as what are the dom- 
inant magnetic effects on low energy cos- 
mic rays. Candidates include magnetic mir- 
roring(discussed in Cesarsky & Volk 1978), 
Alfven weaves(Skilling & Strong 1976), and 
gravitational and turbulence-driven effects. 
Eventually, Fermi acceleration and shock- 
driven acceleration will be added to the model, 
so that the origin and range of these low energy 
cosmic rays can be theoretically explored. 

The main problem that this code addresses 
now is the question of the cosmic ray ion- 



ization rate, and why it has the value that it 
does, connecting it with a flux-spectrum that 
depends on cloud geometry, composition, and 
physical properties like density and electro- 
magnetic properties. At the end of his 2006 
review, Alex Dalgarno stated that "The inter- 
esting question may be not why are [cosmic 
ray ionization rates] so different but why are 
they so similar (Dalgarno 2006)." The prelimi- 
nary results of this study suggest that a combi- 
nation of geometry and magnetic field effects 
may provide the answer to both questions. 

Acknowledgements. We are exceptionally grateful 
to Alexandre Marcowith for organizing the confer- 
ence and surrounding activities, and for his invita- 
tion to the beautiful city of Montpellier. RB.R. also 
thanks the DeMartini Scholarship and Ohio State 
University for its financial support. 

References 

Cesarsky, C. J. & Volk, H. J. 1978, A&A, 70, 
367 

Crank, J., Nicolson, P., & Hartree, D. R. 1947, 

Proceedings of the Cambridge Philosophical 

Society, 43, 50 
Cravens, T. E., Victor, G. A., & Dalgarno, A. 

1975, Planet. Space Sci., 23, 1059 
Dalgarno, A. 2006, Proceedings of the 

National Academy of Science, 103, 12269 
Dalgarno, A., Yan, M., & Liu, W. 1999, ApJS, 

125,237 

Indriolo, N., Geballe, T. R., Oka, T, & McCall, 

B.J. 2007, ApJ, 671, 1736 
Le Petit, P., Nehme, C, Le Bourlot, J., & 

Roueff, E. 2006, ApJS, 164, 506 
McCall, B. J., Huneycutt, A. J., Saykally, R. J., 

et al. 2003, Nature, 422, 500 
Nath, B. B. & Biermann, P L. 1994, MNRAS, 

270, L33 

Padovani, M., Galli, D., & Glassgold, A. E. 

2009, A&A, 501,619 
Rimmer, P., Herbst, E., Morata, O., & RouefF, 

E. 2011, A&A, Submitted 
Skilling, J. 1975, MNRAS, 172, 557 
Skilling, J. & Strong, A. W. 1976, A&A, 53, 

253 

Spitzer, Jr, L. & Tomasko, M. G. 1968, ApJ, 
152, 971 



